Mode resolved travel time statistics for elastic rays in three-dimensional billiards 
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We consider the ray limit of propagating ultrasound waves in three-dimensional bodies made 
from an homogeneous, isotropic, elastic material. Using a Monte Carlo approach, we simulate the 
propagation and proliferation of elastic rays using realistic angle dependent reflection coefficients, 
taking into account mode conversion and ray-splitting. For a few simple geometries, we analyse the 
long time equilibrium distribution focussing on the energy ratio between compressional and shear 
waves. Finally, we study the travel time statistics, i.e. the distribution of the amount of time a 
given trajectory spends as a compressional wave, as compared to the total travel time. These results 
are intimately related to recent elastodynamics experiments on Coda wave interferometry by Lobkis 
and Weaver [Phys. Rev. E 78, 066212 (2008)]. 
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I. INTRODUCTION 

Elastic rays are the fundamental building block of the 
theory of Coda wave interferometry [l[ which over recent 
years has been developed into a well established method 
for the analysis of seismological data @, among oth- 
ers [|| . To the best of our knowledge, only little is known 
about elastic rays as the particle limit of the elastic wave 
equation, even though it is precisely that limit the theory 
in Ref. [l[ relies upon. The main difficulty with elastic 
rays is due to the phenomenon of "ray-splitting" which 
causes the dynamics to become effectively statistical in 
nature Q. 

There is a close analogy in elastic rays and Coda wave 
interferometry on the one hand and the orbits of clas- 
sically chaotic quantum systems and semiclassical the- 
ory in the diagonal approximation [H-Q on the other. 
This analogy is crucial also for measurable quantities 
such as the distortion (introduced in Ref. and scat- 
tering fidelity [§-[l2|. In Ref. Lobkis and Weaver 
started a series of experiments with reverberant ultra- 
sound in three-dimensional Aluminum samples. These 
experiments could be described in terms of the Coda 
wave interferometry (based on the ray picture of a dif- 
fuse wave field), but also in terms of scattering fidelity 
and random matrix theory [Til [l3j. So far the focus 
has been on the form of the fidelity decay but not 
on the overall decay time of the fidelity signal. In or- 
der to explain their results Lobkis and Weaver assumed 
that after a certain transient time, the elastic wave field 
settles on a equilibrium state, where the energy is dis- 
tributed homogeneously and isotropically over the whole 
Aluminum sample, with relative energy shares in shear- 
and compressional waves corresponding to the equiparti- 
tion ratio @, Ell- 
in this paper, we present numerical simulations of the 
propagation of elastic rays in finite three-dimensional 
bodies. Due to ray-splitting an elastic ray spreads into 
an exponentially increasing number of different branches. 



This makes it very difficult to simulate ray-dynamics over 
long times. One of our main achievements consists there- 
fore in the development of an efficient algorithm which 
applies a Monte Carlo approach to deal with these diffi- 
culties. The bodies employed are similar to the ones used 
by Lobkis and Weaver but not identical. Still our sim- 
ulations allow to verify certain assumptions made about 
the wave field and to analyse the effects of eventual vio- 
lations. 

Together with the introduction, this paper is divided 
into 6 sections. The following Sec. [Ill defines our concept 
of an elastic ray. Sec. lIIII describes the random mode con- 
version model introduced in Ref. Q . Our main results are 
described in Sec. IIVI and their relation to experiments is 
discussed in Sec. El We conclude the paper with Scc. lVIl 



II. 



ELASTIC RAYS 



Rays are a well known concept from geometric optics, 
where they are used to approximate the propagation of 
electromagnetic waves in the limit where the wavelength 
is small compared to the typical dimension L of the sys- 
tem. A single ray stands for a transversally concen- 
trated electromagnetic wavepacket which moves through 
an optical system. In the propagation direction, the 
wavepacket may also be localized, but that need not be 
so. Without localization one arrives at a stationary sit- 
uation, where wave energy continuously flows along the 
ray path. 

In the present case, we assume that the wavepacket 
is finite, even in the longitudinal direction. It therefore 
contains a finite amount of energy Q. Since we neglect 
any form of energy dissipation, that amount of energy is 
conserved for all times. However, due to mode conver- 
sion and ray-splitting, the amount of energy Q is shared 
among an increasing number of branches of the elastic 
ray. 

Elastic waves in a homogeneous and isotropic medium 
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come in two different modes: P-waves (compressional 
waves) where the medium undergoes harmonic displace- 
ments in the longitudinal (propagation) direction, and 
S- waves, where the medium undergoes harmonic dis- 
placements in the transversal direction (shear waves). 
Both wave forms have different propagation speeds, Cd 
and c s respectively. For the Aluminum samples studied 
in H [H El: 



Cd = 635 cm/ms 



= 315 cm/ms 



(1) 



As we shall see below, both wave forms convert according 
to specific rules into one-another when the ray is reflected 
at the free surfaces of the body. In the classical ray pic- 
ture this is handled by providing the elastic ray with an 
additional internal degree of freedom, as will be explained 
in detail below, in Sec. Ill Bl 

One should be well aware of the fact that the ray pic- 
ture, when applied to propagating elastic waves in a fi- 
nite solid, breaks down very soon. The transversal size 
of a localized wavepacket increases rapidly in time such 
that after a few reflections, it reaches the extention of 
the whole body. In quantum mechanics this time scale 
would be called the Ehrenfest time JT^]. However, this 
does not mean that the ray picture becomes useless when 
longer times are involved. In the field of quantum chaos 
it is well established that one can construct scmiclassi- 
cal approximations, on nothing else than these classical 
trajectories (Tc| . These approximations may remain valid 
for much longer times - in the case of chaotic two degree- 
of-freedom systems up to times of the order of the Heisen- 
berg time jl9f . Ultimately, the present work might help 
to pave the way towards a similar semiclassical theory 
for elastodynamic systems. 

For the numerical simulations, we use geometries (rect- 
angular block, tetrahedron) with long straight edges, be- 
cause similar bodies have been used in the experiments 
by Lobkis and Weaver and also because of their simplic- 
ity. These bodies however have the disadvantage that 
diffraction on the edges and corners may have consider- 
able effects [13, HJ. The inclusion of diffractive orbits 
will be left to a future work. 



Reflection coefficients 
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FIG. 1. The different reflection coefficients as a function of 
the angle of incidence i9. Rpp (red line), Rsp (green fine), 
Rss (blue line), Rps (pink line). Note, R yx means incident 
i-wave and reflected y-wave. 



Following [22|, when a plane elastic wave hits a free 
surface, we first define a reflection plane. This is the 
plane spanned by the surface normal and the propaga- 
tion direction of the ray. An incident P-wave splits into 
an out-going P-wave and an out-going S'-wave with po- 
larization vector lying in the plane of reflection. The 
reflection amplitudes for both waves are given by: 



Rpp 



Rsp 



sin(2i?) sin(29) - k 2 cos 2 (26) 
sin(2i?) sin(29) + k 2 cos 2 (29) 

2k sin(2?9) cos(29) 
sin(2i?) sin(29) + k 2 cos 2 (29) 



sini? 
sin 9 



(2) 
(3) 



where n = Cd/c s > 1 such that the exit angle 9 is always 
smaller than the entrance angle These angles are al- 
ways measured with respect to the normal of the surface. 
The reflection amplitudes Rpp and Rsp determine the 
amplitudes of the two out-going waves. 

The case of an incident S'-wave is more complicated. 
Before considering the reflection itself, we have to decom- 
pose the wave into one component with polarization in 
the reflection plane (SV-wave) and another with polar- 
ization perpendicular to it (SH-wave). For the SV-wave 
we then have similar reflection coefficients as for the P- 



In general, rays are of a hybrid nature. From a macro- 
scopic perspective, the ray has negligible width and is 
regarded as a one-dimensional object, while, from a mi- 
croscopic perspective, the ray is seen as a plane wave 
which allows to use relatively simple rules to calculate its 
behavior under reflections. Since the system is assumed 
to be macroscopic in size, the shape of the reflecting sur- 
face does not really matter. It is always approximated 
locally as a plane surface. For our purposes it is therefore 
sufficient to consider the reflection laws of plane waves at 
plane surfaces. For simplicity, we shall also assume that 
all reflecting surfaces are free surfaces. 



Rss 



Rps 



sin(20) sin(29') 



3 2 (2tf) 



sin(2i?) sin(29') + k 2 cos 2 (2tf) 

— k sin(4i?) 
sin(2tf) sin(29') + k 2 cos 2 (2i?) 



sin 9' 

sin -d 



(4) 
(5) 



where now 9' > $ such that these equations only apply 
as long as Ksin^ < 1. This introduces the critical angle 
(of incidence) $ cr = arcsin(/t ). However, for SV-waves 
incident at larger angles, the reflected P-wave becomes 
a surface wave with negligible contribution to the wave 
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field, while the reflection amplitude for the SV-wave be- 
comes 

cos 2 (2i?) - 2i/3 sin i? sin(2tf) 
SS ~ ~ cos 2 (2tf) + 2i/3 sin -3 sin(2z?) ' ( ' 

where /3 = y sin 2 i9 — k~ 2 . The absolute value squared of 
i?ss is one in that case, which means that all the energy 
is transferred to the reflecting S-wave. 

SH-waves, i.e. shear waves with polarization direction 
perpendicular to the reflection plane, cannot convert to 
P- waves. They are reflected according to the standard 
law where the reflection angle is equal to the angle of 
incidence. 



B. Classical ray limit 

In order to define the classical ray limit, we consider 
a localized wavepacket of total energy Q. At any mo- 
ment in time this wavepacket has a well defined posi- 
tion r and propagation direction e„. With respect to 
its momentary position r, the wavepacket has a certain 
extension along the propagation direction and perpen- 
dicular to it. Typically, we would imagine a cigar-shaped 
wavepacket, oriented along the propagation direction. In 
the classical limit, we ignore the transversal extension 
of the wavepacket and obtain thereby a one dimensional 
object. For the studies to follow, the wavepacket's ex- 
tension along the propagation direction does not mat- 
ter. Since elastic waves propagate in two different modes, 
the classical description must be extended by some inter- 
nal variables: We need one binary variable to specify 
the mode of the wavepacket which may be longitudinal 
(P) or transversal (S). Moreover, if the wavepacket is 
in transversal mode, we need to record the polarization 
direction by a unit vector e p , perpendicular to the prop- 
agation direction. 

The classical description of elastic waves runs into dif- 
ficulties when it comes to reflections. The problem is ray 
splitting [H, which will be treated as follows: We assume 
that in the vicinity of the wavepacket center, an clastic 
ray may be described as a plane wave. For that plane 
wave, the reflection laws of the previous section [ITX] ap- 
ply. As explained there, the reflected wave will in general 
be split into an 5*-mode and a P-mode branch, propagat- 
ing into different directions. Accordingly, the reflected 
ray splits in two, where the local intensities and polar- 
ization directions in the wavepacket center may be calcu- 
lated from the reflection coefficients of the corresponding 
plane waves. Finally, purely geometric considerations al- 
low to determine the widths of those branches. In that 
way, we obtain the complete information about the two 
reflected branches of the incident ray. Note that an inci- 
dent 5*-mode ray must be considered as a linear combina- 
tion of a SV-wave (polarization direction in the reflection 
plane) and a SP-wavc (polarization direction perpendic- 
ular to the reflection plane). Hence, the P-wave branch 



of the reflected ray is obtained solely from the SV com- 
ponent, wheras the S-wave branch is given as a superpo- 
sition of the reflected S* P-wave and the S-w&ve branch 
of the Sy-component. These relations explain why it is 
necessary to keep track of the polarization direction of 
the S-mode rays. 

For the classical description of the secondary rays, we 
need neither the wavepacket-widths nor the local inten- 
sities, however, we do need the total energy share of each 
branch. In principle, the calculation of the total en- 
ergy share is rather involved since it requires the integra- 
tion of the energy flux across a surface perpendicular to 
the propagation direction, taking into account the vari- 
ation of the wavepacket intensity. Fortunately, for the 
reflected branch without mode conversion, neither the 
wavepacket's intensity profile nor its propagation speed 
change. As a consccucnce, the total energy contained in 
that branch is simply proportional to the intensity cal- 
culated from the corresponding reflection coefficients: 

• An incident P-wave carrying the energy Q will 
transfer the energy R PP Q to the reflected P-modc 
branch. 

• An incident S-wave of the same energy whose po- 
larization direction makes an angle 9 with the re- 
flection plane must be split into the SV-componcnt 
of energy Qy = Q cos 2 9 and the SP-component 
of energy Qu = Q sin 2 9. The total energy of the 
reflected S'-mode branch then is Rg S Qv + Qh- 

Energy conservation now implies that the respective 
mode-converted branch must carry the remaining energy. 
That this is indeed the case, is shown exemplarily for an 
incident P-mode ray in the appendix. 

Returning to the propagation of the clastic ray, we note 
that subsequent reflections lead to further subdivisions 
into an exponentially increasing number of branches. 
Topologically, an elastic ray may be considered as a tree, 
where the initial energy is transported from the trunk to- 
wards the outer (higher order) branches. Even for modest 
travel times, it is soon impossible to keep track of all the 
branches of the elastic ray. 

We therefore employ a Monte Carlo method to sam- 
ple implicitely only over those branches which carry the 
largest amount of energy. For that purpose we trans- 
late energy share into probability. Thus, we replace the 
deterministic energy distribution model by a statistical 
model, where we define an ensemble of rays which all 
start in the same initial state corresponding to a wave 
packet with energy Qq. At each reflection, we choose at 
random, whether the ray follows one branch or the other. 
The probability with which one option or the other is cho- 
sen, are simply given by the relative energy shares cal- 
culated from the corresponding reflection coefficients. At 
the end, the probability to travel through a certain higher 
order branch is given by the product of probabilities for 
the choices made along the history of the given ray. This 
probability agrees with the energy ratio between the to- 
tal energy of the wavepacket in that particular branch 
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and the energy of the initial wavepacket. Thus, in a nu- 
merical simulation, the energy share of a higher order 
branch can be estimated from the number of members in 
the ensemble which terminate in that given branch. 

Our ensemble of rays may as well start in different 
initial states also chosen at random. That is the case, 
when we intend to calculate the evolution of a statisti- 
cal ensemble of initial conditions. In what follows, we 
will consider the following three different types of initial 
conditions: 

(i) Deterministic initial conditions, where we com- 
pletely specify one particular ray, starting at a cer- 
tain point, in a specific direction, either in S- or 
P-modc, and in case with a specific polarization. 

(ii) Surface initial conditions, where we start P-waves 
at a particular point on the surface of the body, 
while the direction of the ray is chosen at random 
within the half sphere pointing into the body. 

(hi) Homogeneous initial conditions, where we start rays 
at random positions inside the body, with random 
directions (isotropic on the whole unit sphere) in 
S- or P-mode with probabilities chosen according 
to the equipartition ratio in Eq. ([7]) and in case 
random polarization direction. 



III. RANDOM MODE CONVERSION MODEL 

In Ref. [|| the authors introduce a simple model, ig- 
noring any geometric effects, where elastic rays undergo 
random and statistically independent mode conversions 
according to the conversion rates a (S-to-P) and (3 (P- 
to-S). Here, the mean free time between two conversions 
is given by a~ l (for the S- mode segments) and (for 
the P-mode segments), respectively. As an additional in- 
gredient, the authors invoke the equipartition principle, 
which states that energy should be distributed equally 
among the different degrees of freedom. Taking into ac- 
count the different wave speeds for P- and S-mode waves, 
the ratio of the energy share between S- and P-mode 
waves is flill 



R = 2 



(7) 



This number also determines the ratio between the con- 
version rates, since the equilibrium condition of the cor- 
responding rate equation demands that on average, the 
number of elastic rays in P-modc and S'-mode are related 
by 



(Ns) 
(N P ) 



= £ = 



(8) 



One last piece of information is necessary in order to 
be able to estimate the values of the conversion rates. 



The authors obtain this from an analogy to room acous- 
tics [lU and a detailed calculation of the mode conversion 
probability during individual reflections The result 



is: 



P = 0.59 c d — . 



(9) 



In order to explain the experiments performed on the 
various Aluminum samples, one requires the statistics of 
tp, the accumulated amount of time a given ray of du- 
ration t spends in P-mode. As an estimate for the mean 
(tp) and its variance var(tp) the authors find 



(tp) 



t 



R+l 



. . 2t 
var(tp) = — 



B 2 



(3 (1 + P)= 



(10) 



The first equation is easily explained on the basis of cr- 
godicity. A more involved calculation is required for the 
second. 

For Fig. [2] we performed numerical simulations for the 
random mode conversion model described above. We use 
random numbers with exponential probability densities 
to generate sequences of alternating P- and S'-mode time 
segments forming a random realisation of an elastic ray. 
With the total time t being fixed we can than measure 
the average P-mode occupation time (tp) as well as its 
variance. These quantities are shown in Fig. [2] as a func- 
tion of t, for three different cases: When the trajectory 
starts with a P-mode segment, with a S-mode segment, 
or by choosing P-mode or S-mode at random acording 
to the equipartition ratio. We find that in all cases, (tp) 
and var (tp) quickly converge to the theoretical values as 
t becomes sufficiently large. 



IV. 



NUMERICAL RESULTS 



In our simulations we use bodies of two different 
shapes, a rectangular block of dimensions 9 cm x 13 cm x 
7.6 cm (giving rise to integrable or possibly pseudo- 
integrable dynamics) and a regular tetrahedron with 
edges of length 10 cm, where the dynamics is probably 
crgodic. In addition, we study each of the two bod- 
ies with and without an internal sphere. That sphere 
is suposcd to provide another free surface for the waves 
(rays) moving inside the body, which renders the dynam- 
ics chaotic. The inner sphere for the rectangle is chosen 
to be relatively large in order to reduce bouncing ball 
orbits. The bodies (including the inner spheres) are de- 
picted in Fig. [3l According to the random mode conver- 
sion model, the volume and the surface area of the bodies 
are important parameters. These are given in table |U 

As explained earlier, the simulation of the classical 
rays is done by launching a large number of rays (up 
to 4 x 10 6 ) with different initial conditions. Those are 
chosen to be either of type (ii) (Surface) or of type (iii) 
(Homogeneous). The former might be considered as cor- 
responding to the experimental situation in Ref. [H, fl3| , 
but if at all this may be true only qualitatively. 
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FIG. 2. Mean and Variance of the P-mode travel time for 
the Lobkis' and Weavers's random conversion model. The P- 
to-S'-conversion rate was chosen as f3 = 162. The different 
colors refer to different initial conditions: starting in S-mode 
(red), starting in P-mode (blue), starting randomly in S- or 
P-mode according to the equipartition ratio (green). In the 
case of the equipartitioned initial condition, the solid green 
line shows the value of the theoretical expectation, Eq. (|10|l . 
Other solid lines (red and blue) show simple rational best fit 
functions to guide the eye. 




FIG. 3. On left a regular tetrahedron with length of 10 cm 
with inner sphere radius 1.2 cm on right a rectangular block 
with length of 9 cm width of 13 cm and height of 7.6 cm with 
inner sphere radius 3.5 cm. 



A. Conversion rates and equipartition ratio 

To verify the accuracy of the random conversion model 
outlined above, we start by analysing the conversion 
rates. For that purpose we perform simulations where the 
initial conditions of the rays are chosen according to the 
expected equilibrium state. Thus, we start rays at ran- 
dom positions inside the body, with random directions, in 
S- or P-mode according to the theoretical equipartition 
ratio (O, and if in S- mode we choose the polarization 
direction also at random. For each ray of pre-defined 
duration t, we then record all periods during which the 
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Volume [cm' 5 ] 


burrace [cm \ 


Rectangular block 


889.2 


568.4 


. . . with inner sphere 


709.6 


722.3 


Regular tetrahedron 


117.8 


173.2 


. . . with inner sphere 


110.6 


191.3 



TABLE I. 
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t/[ms] 

FIG. 4. The conversion rate fi num divided by the travel time 
t, determined from the mean free P-mode occupation time, 
as a function of t, for homogeneous initial conditions, for all 
four different geometries. Tetrahedron with sphere (red line), 
without sphere (green line), Rectangle with sphere (blue line), 
and without sphere (pink line). 



ray happened to travel in P-modc. The average over 
those periods over all rays is just the mean free P-mode 
travel time, and therefore equal to In Fig. 3] the so 
determined conversion rate /3 n um is compared to the the- 
oretical estimate ([9|) for all four geometries considered. 
We find that as soon as t is large enough (for the random 
mode conversion model to become valid, there need to 
occur sufficiently many reflections), /3 num is quite close 
to the theoretical estimate. In fact, as can be observed 
in Fig. [4l /? num never deviates more than 3% from the 
theoretical estimate. Nevertheless, there are systematic 
differences for the different geometries. While (3 nllm ends 
up about 2% above the theoretical estimate in the case 
of the rectangular block with inner sphere, /3, mm ends up 
about 3% below, in the other three cases. We believe 
that this behavior is still acceptable in view of the fact 
that Eq. ([9]) relics on rather rough estimates. 



1. Equipartition ratio 

In the Figs. [S] and O we study the energy share between 
S- and P-mode rays. Since our simulations assume an 
equal amount of energy associated with each ray, that 
energy share can be computed from the relative frequen- 
cies with which we find a particular realization of the 
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FIG. 5. S-mode vs. P-mode energy ratio for samples of dif- 
ferent geometries. Rays are started on the surface in P-mode 
with random directions (Surface type initial conditions). The 
color coding for the different geometries is the same as in 
Fig. [4] The black horizontal line shows the equipartition ra- 
tio R, according to Eq. (J7J). 
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FIG. 6. S-mode vs. P-mode energy ratio for the rectangular 
block with and without internal sphere. The initial conditions 
were chosen homogeneous and isotropic (blue and pink line) 
and starting from a fixed point on the surface (always starting 
in P-mode) (red and green line). The black horizontal line 
shows the value of the equipartition ratio R, Eq. iff}. 



rays in either one of the two possible modes. In Figs. [5] 
and [6l we plot the ratio of these frequencies vs. t, the 
travel time. 

In Fig.[5]the simulation always starts with initial condi- 
tions on the surface [initial conditions of type (ii)], where 
we applied the Monte Carlo sampling with 4 x 10 5 ran- 
dom realizations, for each of the four different geometries. 
Theoretically, we expect that the S/P energy ratio con- 
verges to R = as given in Eq. ([7]). Since in the present 
case, the rays always start in P-mode, the energy ratio at 
small times is close to zero. We observe that two tetra- 
hedrons show a similar behavior, which is very distinct 
from that of the two rectangular blocks. They approach 
the theoretical value for R via damped irregular oscilla- 
tions which are initially very strong. Comparing the two 
tetrahedrons, the presence of the inner sphere, tends to 
reduce the strength of the oscillations such that the the- 
oretical value for R is reached faster. By contrast, the 
rectangular blocks show almost no oscillations at all, and 
quickly go over to a smooth approach of different limit 
values for the energy ratio. In doing so, the rectangle 
with inner sphere comes much closer to the theoretical 
value for R than the rectangle without inner sphere. 

Fig. [6] shows simulations for the two rectangular blocks 
with initial conditions of type (ii) (on the surface) and 
type (iii) (homogeneous). At small times, the energy 
ratio shown starts at zero for initial conditions on the 
surface, and at the value of the equipartition ratio R, 
for homogeneous initial conditions. As t increases, the 
curves first show some minor fluctuations and then start 
to converge to different equilibrium values - in all cases 
clearly below the equipartition ratio. For the rectangular 
block with inner sphere, we observe that independent of 
the type of initial conditions, the energy ratio converges 
to the same equilibrium value. For the rectangle without 



inner sphere, this is not the case, and the equilibrium 
values a very different. 

Extending these simulations up to times of the order of 
5 ms, we have confirmed that the S'/P-cncrgy ratio con- 
verges to finite values in the limit of large times, except 
for the rectangular block without inner sphere. Surpris- 
ingly, we have found that these values do not depend on 
the type of initial conditions applied (Surface or Homo- 
geneous). Concretely, we found the following values: 

{ 16.538 : Tetrahedron without sphere 

16.308 : Tetrahedron with sphere , 

14.685 : Rectangle with sphere ' ^ ' 

9.5/3.2 : Rectangle without sphere 

while the theoretical value is R = 16.384. We can see 
that the rectangular block without inner sphere must be 
considered separately. Only in that case does the equi- 
librium distribution of the energy share depend on the 
initial state (for homogeneous initial conditions we get 
R = 9.5, otherwise R = 3.2). 

B. P-mode occupation time statistics 

In this section, we restrict ourselves to homogeneous 
initial conditions. In our simulations, every trajectory at 
each reflection makes a random choice whether to do a 
mode conversion or not. This results in different trajec- 
tory paths and different amounts of times, the ray spends 
in each mode. In what follows, we study the distribution 
of the P-mode occupation time tp; obviously, the corre- 
sponding S-mode occupation time would provide exactly 
the same information. We have chosen this particular 
quantity because it may be linked to experimental re- 
sults in 0, [l3| with the theoretical models proposed in [l[ 
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FIG. 7. For homogeneous initial conditions, the mean P-mode 
occupation time (tp) for the different geometries, divided by 
the total travel time t. We compare with the theoretical esti- 
mate, Eq. ||T0J, where we have inserted the numerical values 
given in Eq. (|lip . The color coding is as before: Tetrahedron 
with inner sphere (red points), without inner sphere (green 
points), Rectangular block with inner sphere (blue points), 
without inner sphere (pink points). The corresponding es- 
timates for the asymptotic value are depicted as horizontal 
lines of the same color. 



0.3 | r 

0.27 - f 




« ' ' ' ' ' 1 ' 1 ' 1 

0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 

t/M 



FIG. 8. For homogeneous initial conditions, the variance of 
the P-mode occupation times for the different geometries, 
rescaled by /3/t. The color coding for the different geometries 
is the same as in Fig. [7] We compare with the theoretical es- 
timate, Eq. (|1U|I , where /3 has dropped out (black horizontal 
line). Correcting the theoretical expectation for the /3- values 
observed in Fig. 21 leads to the blue horizontal line for the 
rectangular block with inner sphere. 



(Coda wave interferometry) and |8| (random mode con- 
version model). 

We start with the average P-mode occupation time 
(tp), for which the random mode conversion model makes 
the prediction (fTO]) . In Fig. [Ha) we have seen that this 
prediction is indeed very accurate, provided the ensemble 
of rays is in the standard equilibrium state. Fig. shows 
(tp)/t for homogeneous initial conditions which accord- 
ing to theory should converge to (1 + P) _1 at sufficiently 
large times. Indeed for the tetrahedron as well as the 
rectangle with inner sphere, the prediction is fulfilled. 
The ratio (tp)/t converges to the predicted value, if we 
replace the theoretical equipartition ratio with the one 
obtained numerically from our simulations (these values 
are given in Eq. (|11[) ). In the case of the rectangle with- 
out inner sphere (out of range), the convergence of (tp) jt 
is much slower but the limit value still consistent with the 
numerical S/ P energy partitioning. 

Fig. [8] shows the variance of the P-mode occupation 
times, divided by t/f3 with (3 calculated from Eq. ([9]). 
The black solid horizontal line shows the theoretically 
expected value according to the random mode conversion 
model (fTQ|) , where f3 drops out (black horizontal line): 

var(<p) j -> 2 * 2 » 0.1022 . (12) 

Correcting the theoretical expectation for the numeri- 
cally determined /3- values shown in Fig. 3] does not make 
a noticeable difference in the case of the Tetrahedrons 
and leads to the blue horizontal line for the rectangular 



block with inner sphere. In the case of the rectangular 
block without sphere, we find a steep linear increase of 
var(tp) P/t, such that a comparison with a theoretical 
limit value makes no sense. 



In this last figure we find the largest deviations from 
the random mode conversion model. The result for the 
tetrahedron with sphere is about 25% below the theoreti- 
cal value; without sphere, it is about 5% below, while the 
result for the rectangle with sphere is more than 100% 
above the theoretical value. The result for the rectangle 
without sphere is totally off the scale. 

Finally, we show in Fig. [9] the distribution of P-mode 
occupation times tp for trajectories of duration t = 5 ms. 
In the case of the tetrahedrons, we find that the distribu- 
tions are close to Gaussians, with the distribution for the 
tetrahedron with inner sphere being a bit narrower than 
for the tetrahedron without. In the case of the rectan- 
gular blocks, the distributions arc clearly non-Gaussian. 
For the rectangle with inner sphere the distribution is 
strongly asymmetric, while for the rectangle without in- 
ner sphere, the distributions is surprisingly close to a 
Lorentzian. This latter observation can explain the fact 
that the variance of the P-mode occupation times scales 
with t 2 rather than t. since the variance of a Lorentzian 
distribution is infinite. 
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FIG. 9. For homogeneous initial conditions, the histograms 
for the tp statistics for t = 5 ms. The result for the tetrahe- 
drons on the left, with/without inner sphere (red/green solid 
line); for the rectangular blocks on the right with/without 
inner sphere (blue/pink solid lines). 

V. DISCUSSION 
A. Related experimental results 

In Refs. [1,[r|, the authors measured the acoustic long 
time response of short initial ultrasound pulses applied to 
different Aluminum bodies. They studied the cross cor- 
relations between these signals taken at different temper- 
atures. Since the temperature change induces a change 
of the propagation speeds of P- and S'-waves by different 
amounts, a temperature change results in the distortion 
of the acoustic signal and a reduction of the cross correla- 
tions. The reduction of the cross correlations, which has 
also been identified as a scattering fidelity 0, [l(| [H| , can 
be described quantitatively within the theory of Coda 
wave interferometry 

As this theory shows, scattering fidelity (or "distor- 
tion" how that quantity has been called in Ref. Q) is es- 
sentially given by the distribution of P-mode occupation 
times - the quantity studied in the previous Sec. IIVBI 
Assuming Gaussian statistics for these times, the scatter- 
ing fidelity or distortion may be related to the variance 
of the P-mode occupation times as follows: 

D(t) = 2tt 2 (AT ff (S p - 5 S ) 2 var(t p ) , (13) 

where / is the carrier frequency of the elastic wave, AT 
is the change in temperature, while dp = —1.685 x 
10- 4 K-\ and S s = -2.9 x lO^K" 1 arc the thermal 
dilation coefficients for P- waves and S- waves, respec- 
tively [|. 

In Ref. [1], Lobkis and Weaver measured the slope of 
D(t) scaled by AT 2 , / 2 and which according to 

the random mode conversion model should always be the 
same. For our simulations, this is equivalent to deter- 
mining the slope of var(tp) as a function of time, scaled 



by /3 , which should then also give a unique value; see 
Eq. (|12jl . For the Aluminum blocks of different shapes, 
analysed in Ref. 0, the result was a wide spreading of 
experimental values where, the smallest values roughly 
agreed with the theoretical expectation, while the largest 
values where about four times larger [23| . On a qual- 
itative basis, one could observe the tendency that less 
chaotic geometries lead to larger values for the slope. 



B. Our results 

The overall result of the present analysis, shown in 
Fig. [5J is similar in this respect. We also find that the 
slopes can be very different for different geometries. It 
however shows that the relevant quantity is not really 
chaos as measured by Lyapunov exponents but possi- 
bly rather ergodicity. From the four different geometries 
considered, the tetrahedron with inner sphere has the 
smallest slope. It is notable in that case, that the slope 
is about 25% below the theoretical value, which shows 
that the theory doesn't provide a lower bound as one 
might have been conjectured from the experimental re- 
sults. Next comes the tetrahedron without inner sphere. 
Because all its surfaces are plane, the Lyapunov expo- 
nent must be zero in any case. However, the tetrahedron 
has most likely ergodic dynamics. And we find indeed 
that the slope is quite close to the theoretical expecta- 
tion. Only at a considerable distance, we find the rectan- 
gle with inner sphere. The inner sphere clearly leads to 
chaotic dynamics with positive Lyapunov exponents, but 
there are also regions of intcgrablc dynamics and bounc- 
ing ball orbits. That is apparently sufficient to increase 
the slope to more than twice the theoretical value. The 
rectangle without inner sphere is clearly the most reg- 
ular body. For that geometry, the vai(tp) curve rather 
shows a quadratic dependence on t. This is in line with 
the probability density for the P-mode occupation times 
shown in Fig. [HI Its shape is almost a Lorentzian which 
would imply an infinite variance. 

Our simulations have shown that many of the as- 
sumptions within the random mode conversion model are 
rather nicely met. This is true for the P-to-S* conversion 
rate (3 (Fig. |4|), and also for the equipartition ratio if the 
dynamics is ergodic (Fig. [5]). Wc could also confirm that 
the average P-mode occupation time is always related to 
the <S/P-energy ratio (Fig. [7]). Thus, the weak point of 
the random mode conversion model is clearly the in gen- 
eral inaccurate or wrong prediction of the variance of the 
P-mode occupation times. This shows that the succes- 
sion of P-mode and S'-modc segments is in general not 
well described by a Poissonian process. We believe that 
there are two effects coming into play, (i) The duration of 
the individual P-mode segments (these are the shortest 
ones) cannot be really exponentially distributed because 
the rays have to travel a certain distance before having 
the possibility to undergo a mode conversion, (ii) The 
durations of subsequent P-mode and S'-mode segments 
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might be correlated. 

The first effect would lead to a smaller variance of 
the duration of the P-mode segments and thereby to a 
smaller variance of the total P-mode occupation times. 
Thus if the dynamics destroys correlations sufficiently 
rapidly, var(<p) should be smaller than expected theo- 
retically. Our findings suggest that this is indeed the 
case for the tetrahedron with inner sphere. The tetra- 
hedron without inner sphere is ergodic but needs more 
time to destroy correlations. It thus seems that the sec- 
ond effect related to correlations tends to enlarge var(ip). 
For the rectangle with inner sphere, correlations are not 
efficiently destroyed and the variance of the P-mode oc- 
cupation times becomes much larger still. At last, we 
have the rectangle without inner sphere, where var(tp) 
scales with t 2 . 



VI. CONCLUSIONS 

We presented simulations of the propagation of elastic 
waves in three-dimensional bodies of different geometries 
in the limit of classical rays. Because of mode conversion 
at the reflections on the surface of the bodies, one has 
to deal with an exponential proliferation of branches of 
the elastic ray, which is dealt with using Monte Carlo 
sampling. 

Our simulations have shown that there is no unique 
universal equilibrium distribution for elastic rays, at least 
not for bodies with a sufficiently simple geometry. This 
is true even if the dynamics must be considered as com- 
pletely chaotic. For the tetrahedron with internal sphere, 
the equilibrium limit of the S/P energy ratio was close 
but not exactly equal to the theoretical value. Even 
more surprisingly, we found that the homogeneous and 
isotropic distribution of elastic rays with an S/P energy 
ratio equal to the theoretical equipartition ratio need not 
be an equilibrium distribution at all. When chosen as the 
initial condition of an ensemble of elastic rays, we found 
that in the case of the rectangular blocks, the S/P energy 
ratio changes to a different value. 

The main purpose of the present work was to analyse 
some of the conjectures made in Ref. and to investi- 
gate, whether the observed deviations can be explained 
with a model based on classical rays. While many conjec- 
tures and approximations were indeed justified, the most 
problematic assumption was that of random mode con- 
versions with given conversion rates. There, we found 
two counteracting effects: One the one hand the length 
distribution of individual P-mode and S'-modc segments 
is not exponential due to the fact that mode conversions 
can take place only during reflections. On the other hand, 
there are correlations in time where the time scale de- 
pends on the dynamics. In comparison to the random 
mode conversion model, the former tends to reduce the 
variance of the P-mode occupation time while the latter 
leads to an increase. 

Previous publications [TTLfT3j have focused on the form 




FIG. 10. Transformation of the ray widths under mode con- 
version. If the exit angle is equal to the angle of incidence (no 
mode conversion), the width of the reflected ray is unchanged. 
Otherwise the widths are proportional to the cosines of the 
angles ■& and O. 

of the decay function of scattering fidelity which has been 
shown to follow very closely universal random matrix pre- 
dictions. However, our results indicate that the pertur- 
bation strength is an equally interesting quantity. In fact 
it reveals much more system specific information, which 
may even have practical applications, for example in the 
analysis and verification of the properties of mechanical 
components. It would be interesting to design new ex- 
periments, which allowed for a quantitative comparison 
between experiment and numerical simulations. 
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Appendix A: Energy conservation during ray 
splitting 

In Sec. Ill B\ we discuss the splitting of the total energy 
of a reflected wavepacket in the ray limit. In order to de- 
termine the energy share of the mode-converted branch, 
we assumed energy conservation, which is an important 
issue for the self-consistency of our ray model. It implies 
that the reflection coefficients defined in Eqs. (J2j6j) are 
related. Here, we show exemplarily for an incident ray in 
P-mode that the total energy is indeed conserved. 

The amplitude of an elastic ray is physically related to 
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the displacement of infinitesimal volume elements in the 
solid. Therefore, amplitudes up or us (for P-mode and 
S'-mode, respectively) have units of length. The period 
averaged energy density of the respective wave field is 
then given by 



£p,S 



QW 



1 {P,S} 



(Al) 



where g is the density of the medium and w the angular 
frequency of the wave (see Sec. 1.7 of [H[). To demon- 
strate the energy conservation, we show that the energy 
flux is conserved in a quasi-stationary situation where a 
very long wave packet is reflected at a free plane surface. 
It means that the energy flux along the propagation di- 
rection of the incident P-mode ray must be equal to the 
sum of the fluxes along the reflected P-mode and the re- 
flected and mode-converted S-mode branch: In general, 
the amount of energy flowing through a transversal sur- 
face S is given by 



F = / e(f) v n (r) da(r) 
Js 



(A2) 



where v n is the projection of the ray velocity on the sur- 
face normal. If we choose S to be normal to the propa- 
gation direction: 



Fp/s — 



QW 



Cd/ S / u P/s (f) 2 da(f) . (A3) 



Just as the propagation speed, also the energy flux de- 
pends on the wave mode. Now, if the energy flux is really 



conserved, the following relation must hold: 

— c d / u P (r) dcr(r) = — 
J s 



Cd R 



2 

PP 



u'pif) 2 da(f) + c. 



R 2 



SP 



u P (f) 2 da{r) 



(A4) 



where S is a surface perpendicular to the incident ray, 5" 
is a surface perpendicular to the reflected P-mode branch 
and 5"' is a surface perpendicular to the reflected S'-mode 
branch. We have assumed that these surfaces are suffi- 
ciently to the reflection point such that the wave ampli- 
tudes may be considered constant along the line segments 
towards and away from the reflection point. The reflec- 
tion is schematically depicted in Fig.[Tn]scparately for the 
P-mode branch (left hand side) and the S'-mode branch 
(right hand side). Without mode conversion, the geo- 
metrical properties of the ray do not change, such that 
the amplitude u' P relative to S' is just equal to up rela- 
tive to S, which means that the corresponding integrals 
coincide. With mode conversion, the width of the ray 
increases by the factor cos 0/ cos i? > 1. Thus, the inte- 
gral over S" must be scaled by that factor. Therefore, 
dividing by the common factor gw 2 /2, 



c d (1 - R PP ) / up{r) 2 da(r) 
Js 



2 cos Q 

Cs Rsp c^9 



x / u P (f) 2 da(f) (A5) 
Now, since the integrals are the same, we arrive at 



R 2 



pp 



Cd, 



2 cos 8 

Rsp ~^d 



(A6) 



which may be easily verified. Thus the energy flux is 
indeed conserved. 
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